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ABSTRACT 

An  algorithm  is  developed  which  enables  the  computation 
of  horizontal  range  and/ or  depth  from  a  submerged  sound 
source,  using  ray  acoustics  and  the  Lloyd's  mirror 
interference  effect.  The  solution  is  based  on  Snell's  law 
and  involves  integrating  multipath  sound  rays  to  find  the 
difference  in  length  between  the  direct  and  surface 
reflected  sound  paths  from  the  source  to  the  receiver.  This 
difference  in  path  length  is  directly  related  to  the 
observed   Lloyd's    mirror    interference    pattern. 

No  assumptions  as  to  the  mathematical  characteristics  of 
the  sound  velocity  profile  (SVP)  are  made  nor  are  far  field 
approximations  used.  The  solution  is  as  accurate  as  the  SVP 
data   input   to   the    problem. 

A  general  computer  flow  chart  and  basic  language  program 
are  provided  to  allow  local  commands  the  capability  of  pre- 
mission  planning  based  on  specific  operating  area 
environmental    information. 
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I.   INTRODUCTION 

Sound  propagation  between  a  source  and  a  receiver  in  the 
ocean  may  occur  by  many  different  ray  paths,  each  dependent 
on  existing  physical  and  acoustical  conditions.  Crucial  to 
the  determination  of  which  sound  paths  are  of  importance  in 
any  given  problem  are  factors  such  as  the  sound  velocity 
profile,  the  depth  of  the  water,  the  bottom  contour  and 
composition,  the  sea  surface  conditions  and  the  intensity 
and  frequency  of  the  sound. 

This  thesis  deals  with  two  specific  paths  of  sound 
propagation  from  a  source  to  a  receiver.  The  first,  called 
the  direct  path,  does  not  touch  the  surface  or  bottom  of  the 
sea  as  it  travels  from  the  source  to  the  receiver.  The 
second  path  reflects  off  the  sea  surface  before  reaching  the 
receiver  and  is  called  the  reflected  path. 

The  following  explanation  of  the  reception  of  sound  from 

a  single  source  via  these  two  distinct  paths  is  taken 

directly  from  Urick:  [Ref.  1:  pp.  120-123] 

"If  the  sea  surface  were  perfectly  smooth,  it  would 
form  an  almost  perfect  reflector  of  sound.  The  intensity 
of  sound  reflected  from  the  smooth  sea  surface  would  be 
very  nearly  equal  to  that  incident  upon  it." 

He  continues: 

"A  criterion  for  the  roughness  or  smoothness  of  a 
surface  is  given  by  the  Rayleigh  parameter,  defined  as 
R=  kH  sinQ  ; where  k  is  the  wave  number,  (2x)/A  (  A  is  the 
wave   length   of    the   sound    incident   on   the   surface),    H   is 
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the  rms    (sea)    wave  height    (crest  to  trough),   and  6  is  the 
grazing   angle". 


Sea   Surface 


R   =   kH    sin   6 


Figure     1-1. 


Rayleigh    Parameter 
Surface   Roughness. 


(R)     for     Describing     Sea 


"When  R<<1,  the  surface  is  primarily  a  reflector  and 
produces  a  coherent  reflection  at  the  specular  angle  equal 
to  the  angle  of  incidence.  When  R>>1,  the  surface  acts 
as  a  scatterer,  sending  incoherent  energy  in  all 
directions.  When  certain  theoretical  assumptions  are 
made,  (See:  P.  Beckmann,  and  A.  Spizzichino:  Scattering  of 
Electromagnetic  Waves  From  Rough  Surfaces,  p.  93,  The 
Macmillan  Company,  N.Y.,  1963.)  the  (pressure)  amplitude 
reflection  coefficient  n  of  an  irregular  surface  (defined 
as  the  ratio  of  the  reflected  or  coherent  pressure 
amplitude    of    the    return    to    the    incoherent    pressure 


amplitude)     can    be     shown 
expression     m  =    exp    (-R)." 


to    be    given    by    the     simple 


Urick   goes   on   to   explain   the   image    interference. 

"When  the  sea  surface  is  not  too  rough,  it  creates  an 
interference  pattern  in  the  underwater  sound  field.  This 
pattern  is  caused  by  constructive  and  destructive 
interference  between  the  direct  and  surface  reflected 
sound  and  is  called  the  Lloyd  mirror,  or  image 
interference,     effect." 

In   Urick's  sound  propagation  in  the  sea    [Ref.  2:   p.   9-5],  he 

continues   the  discussion   of    surface   reflection   interference: 

"This  sound  field  may  be  divided  into  three  parts,  1)  a 
NEAR  FIELD  close  to  the  source  in  which  the  image  source 
is  too  far  away,  and  the  reflected  sound  is  too  weak  to 
produce  appreciable  interference,  2)  an  INTERFERENCE  FIELD 
in  which  there  are  strong  loops  and  nulls  in  the  signal 
received  by  a  receiver  moving  outward   in   range,   and    3)   a 
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FAR  FIELD  in  which  there  is  an  increasingly  out-of -phase 
condition  between  source  and  image  and  the  intensity  falls 
off  as  the  inverse  fourth  power  of  the  range  (transmission 
loss  increases  as  40  log  r).M 

The    interference     effects    felt    at    the    receiver    is 

expressable    as    the    ratio    of    the     combined    direct    and 

reflected   sound  intensities    to   the  direct  sound    intensity 

alone. 

Stated  mathematically: 

„    ...                Itotal                 I(direct+     reflected) 
Intensity  Ratio   =  =     


or  more  precisely 


■'-direct  -^direct 


2 


2        /RD\  /Rd\ 

Intensity  Ratio   =  1    +     m  |  J  2  y.   I  J  cos(ut) 

\rr/  \rr/ 

Where  RD  is  the  direct  ray  path  length,  RR  is  the  reflected 
path  length,  /x  is  the  surface  pressure  reflection 
coefficient,  u>  (where  u)  =  27rf)  is  the  angular  frequency  of 
the  sound,  and  r  is  the  time  delay  caused  by  the  difference 
in  the  two  path  lengths.  A  formal  derivation  of  the 
intensity  ratio  is  given  in  Appendix  A.  A  more 
comprehensive  treatment  of  the  effects  of  surface  roughness 
on   acoustic   waves    is    given   by   Clay   and   Medwin    [Ref.    3], 
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II.  APPLICATION  of  LLOYD  MIRROR  INTERFERENCE  PHENOMENON  to 
LOCATE  the  HORIZONTAL  RANGE  and  DEPTH  of  the  SOUND 
SOURCE  FROM    the   RECEIVER 

Tactical    tracking   information  from   the   Lloyd's   mirror 

phenomenon  has   been    used   operationally    by   the    Navy   for    many 

years.       In  an   excellent    report   on   the   analysis   of   Lloyd's 

mirror     phenomenon     under     constant     velocity     gradient 

conditions,    Shudde   [Ref.   4:   p.  3],  discusses  the  difficulty 

of    solving    the    surface    reflection    interference    problem. 

Accurate   error   estimation   for   any   acoustical  condition   other 

than   the  isovelocity   case   is    not   practical   unless   a    large 

computer     is     used.         Since     Shudde's     1977     report,      the 

advancements  in  memory  size,    and   speed  of   operation  of   small 

portable    desk    top    computers    now    allow    an    accurate    solution 

of  source  position  for  all  sound  velocity  conditions.      The 

following    mathematical    solution    is    based    on    the    lengths    of 

the  sound  paths  determined  from   Snell's  Law.      As  explained 

by  Urick,    [Ref.   1:   p.   116], 

"One  of  the  most  important  practical  results  of  ray 
theory  is  Snell's  Law,  which  describes  the  refraction  of 
sound   rays    in   a  medium    of   variable   velocity." 

The    solution   is   not  dependent   on  any   isovelocity   or   constant 

gradient    conditions    other    than   the   reasonable   assumption 

that  no   horizontal    sound  speed   gradients  exist   throughout 

the    operating    area.       The    mathematical    solution   presented    in 

Section    III    makes    no    assumption    as    to    the    mathematical 
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characteristics  of  the  sound  velocity  profile  nor  does  it 
use  far  field  approximations.  Past  solutions  depended  on 
specific  sound  velocity  profiles,  such  as  an  isovelocity  or 
"straight  line"  profile  which  resulted  in  significant 
positional  errors  for  some  tactical  applications  when  the 
actual  sound  velocity  profile  did  not  match  the  predicted 
model . 

In  applying  the  solution  to  the  Lloyd's  mirror 
interference  effect  one  must  understand  the  relationship  of 
the  propagating  sound  to  the  interference  pattern.  Coppens, 
Sanders  and  Dahl    [Ref.   5    :  p.  55  ]   explained, 

"...the  noise  radiated  by  underwater  vehicles  consists 
of  two  contributions:  broadband  noise  (similar  to  the 
hissing  found  between  FM  radio  stations)  and  tonals 
(single  frequency  tones,  like  whistles).  Thus,  the 
radiated  noise  consists  of  all  frequencies  between  certain 
broad  limits  with  certain  specific  frequencies  occurring 
with  high  intensity.  Each  of  these  frequencies,  whether 
of  broadband  noise  or  tonal,  has  its  own  interference 
pattern. 

Because  each  frequency  has  its  own  period  and 
wavelength,  each  will  have  a  different  phase  delay  for 
propagation  over  the  direct  and  surface  reflected  paths. 
As  a  result,  the  ranges  between  source  and  receiver  for 
which  the  surface  reflected  signal  tends  to  cancel  the 
direct  signal  will  be  different  for  each  different 
frequency.  Thus,  the  ranges  at  which  different 
frequencies  will  have  nulls  at  the  receiver  will  be 
different.  This  means  that  as  the  range  between  source 
and  receiver  changes,  the  frequencies  which  have  nulls  at 
the   receiver    will   also    change." 

Figure  II-1    shows  a    schematic  time  vs.    frequency   plot   of    the 

interference    pattern    as    the    distance   between    the    source    and 

receiver    is    decreased    and    then    increased.       The    distinct 

tonals    will    most    readily    show    the    effects    of    Doppler    shift 
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(apparent  frequency  change)  as  minimum  separation  range,  is 
passed  (see  t  =  0  in  Figure  II-1).  The  light  (unshaded) 
hyperbolic  shaped  regions  centered  about  t  =  0  are  known  as 
the  Lloyd's  mirror  interference  pattern  for  the  broad  band 
noise.  The  frequency  values  F-j  ,  ?2i  F3>  F4  represent  the 
minimum  frequency  of  each  hyperbola.  The  actual  frequency 
from  the  sound  source  that  is  being  nullified  by  the  surface 
reflection  interface  varies  with  time  as  the  hyperbolas  are 
traced  out.  The  measured  separation  in  hertz  between  any 
two  successive  nulls  at  any  time  t,  is  termed  the  difference 
frequency,  A  f ,  and  is  the  frequency  value  required  in  the 
problem  solution  discussed  in  the  following  Section.  (See 
Section  III.D.2) 
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Figure    II-1.      Frequency   vs.    Time    Display  of  Lloyd's    Mirror 

Interference   Pattern 
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III.   RANGE  vs.  DEPTH  SOLUTION  USING  the  LLOYD ' S  MIRROR 

INTERFERENCE   PHENOMENON 

A.   THEORETICAL  BACKGROUND 

1 .   Snell ' s  Law 

Snell's  Law,  [Ref.  1:  p.  116,  Ref.  6:  p.  120],  is  a 
result  of  Fermat's  minimum  time  principle  which  states  that 
each  ray  will  travel  the  path  requiring  minimum  time  of 
travel    from   source   to   receiver. 

Mathematically,  Snell's  Law  is  expressed  as: 


Snell  Constant  = 


sin  0 
c(z) 


© 


Where  0  is  the  angle  that  the  ray  makes  with  the  vertical 
and  c(z)  is  the  sound  velocity  for  a  given  depth  z.  (See 
Figure    III-1  ) 

Rinitial 
Surface 


zinitial~ 

Depth  ( z  ) , 

Meters     + 


origin 


+  Range , 
Meters 


Source 


Raypath  direction  at  Depth  Z^, 

Range  R^ 


Figure  III-1.      Snell's  Law  Applied  to    Initial  Sound  Source 

Ray  Angle. 
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The  sound  ray  maintains  the  same  Snell  Constant  value 
throughout  its  transmission  path.  Sound  velocity  changes 
result   in    an   increasing    or   decreasing   angle   0. 

2.      Ray   Path   Integration 

Expanding  the  diagram  of  the  sound  ray,  the 
differential  raypath  length,  ds,  can  be  described  as  a 
function  of  the  differential  change  in  depth,  dz,  and  the 
differential   change   in  horizontal   range,    dr.       (Figure   III-2) 


Sound 
Ray 


Source 


ds    =      \/dz2    +   dr2 


Figure     III-2.        Pythagorean    Theorem    Applied    to    Ray    Path 

Length. 


To  calculate  the  actual  path  length  of  the  sound 
ray,  smaller  and  smaller  segments  of  path  length,  As, 
horizontal  range,  Ar,  and  depth,  Az  are  used  until  the 
differential  segments  ds,  dr,  and  dz  are  approached.  This 
tends  to  smooth  the  curve  of  the  ray  into  a  sequence  of 
straight  line  segments.  The  path  length  traveled  by  a 
particular  ray  can  then  be  determined  by  integrating  between 
the  source  depth  and  receiver  depth.  The  possibility  of  the 
direct    ray   crossing   the    same    depth    more    than   once    is    not 
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considered  a  problem  in  this  solution  due  to  the  short 
ranges  involved  in  Lloyd's  mirror  observations  and  the  long 
ranges  required  for  a  sound  ray  to  curve  back  through  the 
same  depth  before  reflecting  off  the  surface.  To  handle  the 
fact  that  the  surface  reflected  ray  does  cross  the  same 
depth  twice  in  the  ranges  pertinent  to  this  discussion  and 
to  avoid  errors  in  integration,  the  concept  of  a  geometric 
equivalent  image  [Ref.  1:  p.  123],  of  the  source  is  used. 
The  image  reflected  path  and  the  actual  reflected  path  are 
the   same    length.    (Figure    III-3) 


Image  Reflected  Path 


c(z) 
SVP 


Surface      o 


+  Range 
Meters 


Actual 

Reflected 

Path 


+  Depth  * 
Meters 


Figure    III-3.      Actual   and  Image   Sound  Paths   From   Source  to 

Receiver  for  a    Given  Sound  Velocity   Profile 
( SVP ) . 


The  curved  sound  path  shown  in  Figures  III -3  and 
III-4  is  due  to  the  change  (in  general)  in  the  sound 
velocity  as  the  depth  changes.  It  is  possible  to  show  that 
the  path  from  A   to   B  is  the    same    length   as    the   path   from    A 
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to  B  .  Also,  by  inspection  of  the  sketch,  the  horizontal 
range  from  A  to  B  or  A  to  B  is  the  same.  Therefore,  for 
the  reflected  path  calculations  of  path  length  and 
horizontal  range,  the  image  reflection  depth  of  -z^  will  be 
used. 

B.       HORIZONTAL    RANGE    SOLUTION 
1 .      Direct   Path 

The  horizontal  range  traversed  by  the  direct  path 
sound  ray  can  be  calculated  from  Snell's  Law  and  the 
Pythagorean  Theorem  for  right  triangles  applied  to  Figure 
III-2. 


From    (?)     ,    the   Snell   Constant,    K1    is 

sin  0^ 
c(zi) 


K1    =  .  (1A) 


This  is  the  Snell  constant  for  the  particular  direct  path 

ray  that  ultimately  travels  from  the  source  to  the  receiver, 

but  initially  starts  at  depth  z^  and  moves  in  direction  <>^. 

Using  the  Pythagorean  Theorem  on  Figure  III-2 

dr 
sin  0  =      ____ (Z) 


(dr)2  +  (dz)2 


Solving         HA)       for   sin  0,        where  0   =  variable  for   the 
^■"^  particular    direct    path    ray 

tagged    with    Snell    constant 

K1 '  /~N 

sin  0   =   K1    c(z)  (3) 
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setting  (J)  and  (3)  equal  to  each  other  yields: 


dr 
K1  c(z)  =   (?) 


(dr)2  +  (dz)2 


The  right  side  of  (7)  can  be  simplified  by  dividing 
numerator  and  denominator  by  dr. 


K-,  c(z)  = 


'  *  (S) 


Solving  for  dz/dr  follows  in  steps 


1  + 


1   + 


0 


K-,  c(z) 


[K-,  c(z)]2 


/dz\  2  1    -    [K1    c(z)]2 

\dr/  [K1    c(z)]2 


dz  /1    -    [K-,    c(z)]2 

dr  W    [K-,    c(z) ]2 
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=      + 


dz  /1    -    [K-    c(z)]2 

© 


dr  -  K-,    c(z) 

Now,     solving    for    dr,     the    incremental    horizontal 
range   yields: 

^    c(z) 


dr  =      + 


y/l    -    [K-,    c(z)]2 


dz  0 


Integrating     0     with    respect    to    depth,     with    the 

limits   of    integration  being    the   depth    of    the    source   and    the 

depth   of    the   receiver,     yields: 

Source 
Depth 


K«!    c(z) 
R    (direct   path)    =      +   I    dz  (7) 


VI    -    [K1    c(z)]2 

Actual 

Receiver 

Depth 

(along   path   of    ray 

"tagged"   K1 ) 

The  "Range"  is  always  positive  and  thus  the 
appropriate  +  or  -  sign  is  used  to  make  this  integral  always 
positive  regardless  of  the  integration  limits.  For  ease  of 
discussion  in  the  rest  of  this  solution,  the  receiver  will 
be  assumed   to   be  at   a   shallower   depth. 
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2.   Reflected  Path 

Recalling  the  image  depth  discussion  and  noting  a 

different  angle,  6^,  for  the  reflected  ray  is  required 

(Figure  III-4),  which  results  in  a  different  Snell  constant, 

K2,  where 

sin  {9^) 

K2  =   .  (8) 

c(Zi) 

It  follows  directly  that  (r)     can  be  rewritten  for  the 
reflected  horizontal  range  as: 


Source 
Depth 

K2  c(z) 
R  (reflected  path)  =1 dz  (9) 


\/l  -  [K2  c(z)]2 


Image 

Receiver 

Depth 
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-Zi 


+  z. 


Depth 
Meters 


+  ♦ 


R(  reflected  path)- 

Image  Receiver 


+  Range 
••  Meters 


c(za) 


R( direct  path) 


Figure  III-4.  Common  Direct  and  Reflected  Sound  Paths  From 

Source  to  Receiver. 


C.   SOUND  RAY  PATH  LENGTH  SOLUTION 
1 .   Direct  Path 

The  direct  path  ray  length  can  be  calculated  by 
returning  to  Figure  III-2,  applying  the  Pythagorean  Theorem 
and  noting: 


ds   = 


(dr)2  +  (dz)2 


realizing  the  path  length  S  from  receiver  to  source  is  given 
by: 


Source 
Depth 


Source 
Depth 


S  = 


/ 


ds 


1  V  (dr)2  + 


(dz) 


Receiver 
Depth 


Receiver 
Depth 
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dz 
Which  can  be  simplified  by    multiplying    by        


\l  dz2 

as   follows: 

Source 
Depth 


S    (Direct   Path)    =  /      /1    +     /dr\   2 

Receiver 
Depth 


/dr\ 

Substituting  the  equivalent  of  I  —  I  from  the  reciprocal  of 

©, 

Source 
Depth 


(Direct  Path)  =  K-, 2  c(z)2 


1  +  - dz 

1  -  K^2   c(z)2 


Receiver 
Depth 


Source 
Depth 


S  (Direct  Path)  =   f  /1    K12c(z)2+K12c(z)2 

dz 


1  -  [*.,  c(z)]2 


Receiver 
Depth 
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Source 
Depth 


S  (Direct  Path)  =    I   dz        (12 


/ 


1  -  [K,  c(z)] 


2 


Actual 

Receiver 

Depth 


2.   Reflected  Path 


Using  the  same  arguments  about  the  reflected  path  as 
presented  for  horizontal  range  and  substituting  K2  for  Ki 
and  image  receiver  depth  into  Ml)  yields: 

Source 

Depth 

dz        (Q) 

\Jl    -  [K2  c(z)]2 


Image 

Receiver 

Depth 


D.   EXPLANATION  OF  SOLUTION 
1 .   Some  Observations 

Equations  (T)  and  ui)  are  path  integrals  for  the 
direct  ray,  using  the  same  Snell  constant  and  integration 
limits.  Further,  equations  (s)  and  (V3)  are  path  integrals 
of  the  reflected  ray  and  thus  use  a  different  Snell  constant 
from  the  direct  path  integrals,  but  use  the  same  integration 
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limits   respectively   for  the   horizontal   range  and   ray  path 
lengths. 

2.      Solution 

The  key  to  the  eventual  numerical  solution  is  the 
equality  between  the  direct  path  horizontal  range  and  the 
reflected  path  horizontal  range.  Although  the  actual  path 
length  of  the  reflected  ray  is  longer  than  the  direct  path 
ray  between  common  source  and  receiver,  the  horizontal 
distance   traveled  by  both   is   the    same. 

By  calculating  all  of  the  initial  ray  angles  that 
are  possible  for  each  depth  of  the  source  of  interest,  two 
horizontal  range  arrays,  one  for  the  direct  path  and  one  for 
the  reflected  path,  are  formed.  Each  of  these  arrays  is 
linked  respectively  to  a  direct  and  reflected  ray  path 
length  array  by  the  corresponding  Snell  Constant  used  to 
determine  both  the  horizontal  range  and  sound  ray  length  for 
each  path,  (See  Figure  III-4).  For  example,  suppose  the 
initial  direct  path  ray  angle  of  15  degrees  from  the 
vertical  covered  a  horizontal  range  of  20  meters  in 
traveling  from  the  source  to  the  receiver.  Similarly,  the 
initial  reflected  path  ray  angle  of  8  degrees  from  the 
vertical  covered  the  same  horizontal  range  of  20  meters  from 
the  source  to  the  receiver.  Then  the  difference  between  the 
reflected  ray  path  length  for  8  degrees  and  the  direct  ray 
path   length   for    15   degrees    could  be    computed    (Figure    III-5). 
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a  • 


+ Depth 
Meters  * 


Receiver, 


Range 
Meters 


8  degrees 
1 5  degrees 

Source 


Figure  III-5.   Example  of  Reflected  Path  Ray  Interfering 

With  Direct  Path  Ray  at  the  Receiver. 


It  is  this  difference  in  ray  path  length, 
AS  =  S  -  S  ,  for  a  wide  band  of  frequencies  that  causes  the 
time  delayed  phase  shift  sensed  at  the  receiver  and  known  as 
the  Lloyd's  mirror  interference  phenomenon.  Since  both  the 
direct  and  reflected  rays  are  traveling  at  the  same  sound 
velocity  at  the  receiver,  the  ray  path  length  difference  can 
be  directly  converted  to  the  frequency  difference,  Af  as 
observed  at  the  receiver  and  seen  in  the  Lloyd's  mirror 
interference  pattern. 

The  intensity  ratio  at  the  receiver  is  given  by 


Intensity  ratio  =   1  +  m 


©2  -  > id 


COS  (  a;  t  ) 
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The  ratio  reaches  a  minimum  value  corresponding  to 

maximum  interference  when   cos(u  r)  =1.   This   occurs   when 

w  t  =  2  t  or  any  integer  multiple  of  2  ir   . 

Note  that    X  =   c/f   where  X  =  wavelength  (m) 

c  =  phase  speed  (m/s) 
f  =  frequency  in  (hz) 

k  =  wave  number  (m   ) 

and  k  =  2  *■  /  X 

Now  let  A0  equal  the  phase  difference  seen  at  the 
receiver  due  to  the  two  path  length  differences  and  defined 
as, 

A0  =  k  (SR-  SD),  where 

SR  =  Path  length  of  the  reflected  ray  and 

SD  =  Path  length  of  the  direct  ray. 

Let  A  S  =  (SR-  SD) 

So  that  the  first  null  occurs  when 

A0  =  k(AS)=2T=u>T 

Sok  (AS)  =wt=  2r=  2  *■  (  A  S )  /  X  or 
A  S  =  2tX/2t   =  X 
therefore  AS  =  c/f-j 

so  the  first  null  occurs  at  frequency 
f<l    =   c/  AS 

The  second  null  occurs  at 

w  t  =  4  r   so  that 

AS    =    47rX/27r=      2X=      2c/f2 
therefore    second  null   occurs   at   frequency 
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f2  =  2c/A  S   =  2(f1  ) 
which  shows  that  f2  -  f-j  =  f-j 
and  so  on.   So  that  for  the  nth  null 

fn  =  nc/d  S   =  nf"i 
and      fn  -  fn_i  =  f-j  or  more  generally 

fm   ~   fn  c 


=   f1    = 


m   -   n  AS 

where  m  >>  n  for  best  accuracy  so  that  the  measured 
difference  in  frequency  between  any  two  successive  nulls  can 
be  related  to  the  difference  in  path  length  required  to 
produce   that   null   by   the   equation 


C( receiver) 


AS 


where   c(receiver)    is    the   actual    sound   velocity   at 
the   receiver   at    the   time   of   reception. 
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IV.      COMPARISON   of    the   ANALYTICAL   and   COMPUTER    SOLUTIONS    for 
the    ISOVELOCITY    SOUND   PROPAGATION    CASE 

This  analysis  demonstrates  the  accuracy  of  the  computer 
solution  by  using  an  isovelocity  sound  propagation  model. 
The  constant  sound  velocity  condition  results  in  straight 
line  transmission  paths  as  the  sound  travels  from  the  source 
to  the  receiver.  By  solving  the  geometric  model  analytically 
and  comparing  the  results  to  the  computer  output  one  can  see 
the  accuracy  of  the  Romberg  integration  technique  coupled 
with  the  four  integrals  derived  in  section  III  .  The 
isovelocity  case  is  the  simplest  propagation  model  yet  it 
exercises  all  phases  of  the  computer  program's  main 
calculation  algorithm.  The  only  section  of  the  overall 
program  not  tested  by  the  isovelocity  case  are  the  sound 
velocity  calculations  performed  as  the  user  inputs  a  varying 
temperature  profile.  Table  IV- 1  compares  the  difference  in 
path  lengths  between  the  direct  and  reflected  sound  paths  as 
calculated  analytically  and  by  the  computer.  Of  the  20  data 
points  compared,  the  computer  missed  the  exact  ray  path 
length  for  a  horizontal  range  of  750  meters  by  only  19.16 
centimeters  in  the  worst  case  and  averaged  5.02  centimeters 
error  over  all  20  cases.  Table  IV-2  shows  the  resulting 
horizontal    range   error  due   to   the   error   of   the   computer's 
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calculation  presented  in  Table  IV-1.  The  average  range  error 
for   the    20    points    is    only   1.52    meters. 

The  computer  program  steps  through  sea  depths  of  0  to 
200  meters  in  10  meter  increments  where  each  depth 
represents  a  possible  target  depth.  To  calculate  a  complete 
set  of  Snell  constants,  initial  sound  ray  angles  of  3  thru 
177  degrees  in  increments  of  3  degrees  (  measured  from  the 
vertical  )  are  used.  Ultimately  the  integration  section 
derives  a  path  length  difference  between  the  surface 
reflected  ray  and  the  direct  path  ray  for  each  of  these 
selected  starting  ray  directions. The  59  initial  rays  times 
the  20  depth  increments  yields  a  total  of  1,1 80  data  points. 
These  point  solutions  are  used  to  interpolate  values  of  path 
length  difference  for  points  every  10  meters  of  horizontal 
range  out  to  1000  meters  which  results  in  a  two  dimensional 
matrix  containing  2000  (  20  depth  increments  times  100  range 
increments    )    horizontal   range  vs   depth   elements. 

To  compare  the  computer  solution  to  actual  geometric 
results,  five  ranges  for  each  of  four  depths  were  chosen  at 
random.  In  all  comparisons  the  receiver  is  set  at  a  depth  of 
27    meters    (approximately    90    feet). 
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TABLE    IV- 1 


SOURCE 


ISOVELOCTTY    CASE 
COMPARISON    BETWEEN 
GEOMETRIC  AND     COMPUTER    SOLUTIONS 


PATH    LENGTH    DIFFERENCE 

ANALYTICAL 

COMPUTER 


DEPTHS 
METERS 

HORIZONTAL   RANGES, 

METERS 

20 

100 

380 

750 

10 

15.8107 

5.1908 

1.4170 

0.7195 

15.8228 

5.1939 

1.4418 

.6693 

50 

49.0755 

23.5992 

7.0274 

3.5897 

49.0196 

23.6220 

7.0389 

3.5714 

90 

52.5987 

35.7218 

12.4171 

6.4289 

52.6316 

35.7143 

12.4688 

6.2267 

130 

53.3450 

42.5840 

17.4439 

9.2169 

53.3808 

42.6136 

17.5644 

9.0253 

170 

53.6208 

46.4312 

22.0132 

11 .9302 

53.5714 

46.4396 

22.0264 

11 .7739 
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TABLE  IV-2 
COMPARISON  of  PROGRAM  RANGES  vs.  ACTUAL  RANGES 


SOURCE 

DEPTHS 
METERS 


ACTUAL 


HORIZONTAL   RANGES,  METERS 


20 


100 


380 


750 


10 

50 
90 
130 
170 


19.954      99.935     373.42 


20.134      99.875     379.37 


19.42 


21  .26 


806.3 


753. 87 


19.745  100.038  378.33  774.84 


99.815  377.075  766.41 


99.93  379.725  760.475 


COMPUTER 
HORIZONTAL      RANGES,    METERS 
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V.       COMPUTER    PROGRAM   FLOW    CHART 

A.      FLOW    CHART   and   PROGRAM    BACKGROUND 

The  following  block  diagram  flow  chart  demonstrates  the 
specific  steps  required  to  obtain  an  accurate  horizontal 
range  vs.  depth  solution  using  the  Lloyd  mirror  interference 
phenomenon.  Specific  values  or  data  used  in  the  included 
program  will  be  given  at  the  bottom  of  each  appropriate 
block. 

The  actual  numerical  integration  section  of  the  program 
is  taken  directly  from  the  Romberg  integration  program 
provided   by  Miller    [Ref.    7:    pp.    215-218]. 

The  sound  velocity  equation  used  in  the  program  was 
developed  by  Coppens  [Refs.  6  and  8]  and  selected  after 
several  equations  of  varying  complexity  were  considered. 
Specifically,  equations  derived  by  Wilson  and  Leroy  [Ref.  1: 
p.  105],  plus  Mackenzie  [Ref.  9]  were  compared  with  Coppens. 
All  results  were  within  1  to  2  centimeters  per  second  of 
each  other.  Coppen's  equation  was  chosen  because  of  the 
ability  to  correct  for  not  only  temperature,  pressure  and 
salinity,    but    also,    latitude    variations. 
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B. 


FLOW  CHART 

1 .   Initialization  and  Main  Menu 


1. 
2. 
3. 

INITIALIZATION 
DEFINE  ALL  VARIABLES 
DIMENSION  ALL  ARRAYS 
DEFINE  TOLERANCE 
DESIRED  FOR  NUMERICAL 
INTEGRATION 

" 

RANDOM  ACCESS  FILE 

CREATE  A  RANDOM 

ACCESS  FILE  TO 

STORE  THE  RANGE 

VS.  DEPTH  MATRIX 

SOLUTION 

t. 

? 

RETURN  FROM 
ERROR  STATEMENT 
FOLLOWING 
INTEGRATION  ERROR 

RETURN  TO 
MENU  FROM 
PROGRAM 
SUBROUTINES 
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DISPLAY  MENU  OF 

PROGRAM  SUBROUTINES 

ALLOWING  USER  TO 

CHOOSE  DESIRED 

OPERATION  AS 

FOLLOWS: 


1 .   REINITIALIZE  RANDOM 
ACCESS  ANSWER  FILE 
AFTER  INITIAL 
RUN  AND  PRIOR  TO 
ALL  SUBSEQUENT  RUNS 


INPUT  WATER  TEMPERATURE 
INFO  AS  A  FUNCTION  OF 
DEPTH 


3.   INPUT  SEA  SURFACE  TEMP- 
ERATURE AND  AVERAGE 
GRADIENT 


4.  INPUT  SOUND  VELOCITY 
PROFILE  AS  A  FUNCTION 
OF  DEPTH 

5.  COMPUTE  HORIZONTAL  RANGE 
VS.  DEPTH  SOLUTION  AND 
STORE  IN  ANSWER  FILE 


6.  DISPLAY  ANSWER  FILE  ON 
THE  CRT  SCREEN  OR  LIST 
ON  A  PRINTER 


7.   EXIT  PROGRAM  AND  RETURN 
TO  OPERATING  SYSTEM 


<D 


S2) 


\& 


<!) 


-0 
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2.   Reinitialize  Random  Access  Answer  File 


G> 


WARN  USER  THAT 
THIS  SECTION  WILL 
DESTROY  CURRENT 
ANSWER  FILE 

QUERY  USER 
REINITIALIZATION 
DESIRED,  YES  OR  NO 


NO 


RETURN 
TO  MENU 


® 


YES 


1  . 

~> 
—  • 

3. 


REINITIALIZE  FILE 

CLOSE  ANSWER  FILE 

KILL  (DELETE)  ANSWER  FILE 

OPEN  NEW  ANSWER  FILE 


RETURN 
TO  MENU 
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3.   Input  Water  Temperature. Profile 


INPUT  SALINITY  IN  PARTS 
PER  THOUSAND  FOR  OPERATING 
AREA,  ENTER  A  VALUE  OF  35 

IF  SALINITY  IS  UNKNOWN 


INPUT  LATITUDE  IN  WHOLE 
DEGREES  FOR  OPERATING  AREA 


INPUT  WATER  TEMPERATURE  IN 
DEGREES  CELSIUS  AT  UNIFORM 
INCREMENTS  OF  DEPTH  IN  METERS 
BEGINNING  WITH  SURFACE  (PROGRAM 
STEPS  IN  5  METER  INCREMENTS 
DOWN  TO  200  METERS) 
CALCULATE  SOUND  VELOCITY  IN  METERS 
PER  SECOND  FOR  EACH  1  METER 
INCREMENT  OF  DEPTH  AND  STORE 
ANSWER  IN  TEMPORARY  ARRAY 


RETURN  TO 
MENU 
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4.   Input  Sea  Surface  Temperature  And  Average  Gradient 


© 


INPUT  SALINITY  IN 
PARTS  PER  THOUSAND  FOR 
OPERATING  AREA,  ENTER  A 
VALUE  OF  35  IF  SALINITY 
IS  UNKNOWN 


r— ■ 

INPUT  LATITUDE  IN  WHOLE 

DEGREES  FOP  OPERATING  AREA 


T 


INPUT  SEA  SURFACE 
TEMPERATURE  IN 
DEGREES  CELSIUS 


INPUT  AVERAGE  TEMPERATURE 
GRADIENT  IN  DEGREES 
CELSIUS  PER  METER  DOWN 
|TO  LOWEST  DEPTH  OF  INTEREST 
(PROGRAM  USES  200  METERS 
AS  LOWEST  DEPTH ) 


CALCULATE  SOUND  VELOCITY 
IN  METERS  PER  SECOND  FOR 
EACH  1  METER  INCREMENT  OF 
DEPTH  AND  STORE  IN 
TEMPORARY  ARRAY   


RETURN 
TO  MENU 


<J> 
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5.   Input  Sound  Velocity  Profile 


© 


INPUT  SOUND  VELOCITY 
IN  METERS  PER  SECOND  FOR 
INCREMENTS  OF  DEPTH  IN  METERS 
(PROGRAM  USED  5  METER 
INCREMENTS)  STORE  VELOCITIES 
IN  A  TEMPORARY  ARRAY 
(PROGRAM  INTERPOLATES  1  METER 
INCREMENTS  FROM  5  METER  INPUTS) 


j  RETURN 
I  TO  MENU 


0 
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6.   Compute  Range  vs.  Depth  Solution  And  Store  In  Random 
Access  File 


<3> 


INPUT  HYDROPHONE/RECEIVER 
DEPTH  IN  EITHER 
FEET  OR  METERS 


|  LOOP  THROUGH  EACH  POSSIBLE 
! SOURCE  DEPTH,  (THIS  PROGRAM 
I  USES  1 0  METER  INCREMENTS  OF 
DEPTH  DOWN  TO  200  METERS) 
COMPUTING  PATH  LENGTH 
DIFFERENCES  FOR  EACH 
INCREMENT  OF  HORIZONTAL 
RANGE,  (THIS  PROGRAM  USES 
10  METER  INCREMENTS  OF 
HORIZONTAL  RANGE  OUT 
TO  1000  METERS) 
BY  DOING  THE  FOLLOWING: 

1  .   DETERMINE  THE  UPPER 

AND  LOWER  LIMITS  FOR  THE 
FOUR  SEPARATE  NUMERICAL 
INTEGRATIONS 

2.   LOOP  THROUGH  EACH 

POSSIBLE  INITIAL  RAY 

ANGLE  FROM  EACH  SOURCE 

DEPTH  COMPUTING  THE  FOLLOWING: 

a.  RETRIEVE  THE 
REQUIRED  SOUND  VELOCITIES 
FOR  EACH  INTEGRATION 

b.  COMPUTE  EACH 
UNIQUE  SNELL  CONSTANT 
GIVEN  AN  INITIAL  RAY 
ANGLE  AND  SOUND  VELOCITY 
FOR  THE  SOURCE  DEPTH 

C.   NUMERICALLY  SOLVE 

THE  FOUR  INTEGRALS  DERIVED 
IN  THE  MATHEMATICAL 
SOLUTION  DISCUSSION, 

tcont.  p.  9>- — 
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-(cont.  from  p.  8) 

USING  A  ROMBERG  IMPROVE- 
MENT [REF.  7]  ON  THE 
TRAPAZOIDAL  NUMERICAL 
INTEGRATION  TECHNIQUE 

d.  IF  INTEGRATION  ERROR 
OCCURS  WARN  USER  AND 
RETURN  TO  MENU 

e.  STORE  EACH  INTEGRATION 
ANSWER  IN  THE  APPROPRIATE 
TEMPORARY  ARRAY 

f.  CONTINUE  TO  LOOP  THROUGH 
ALL  POSSIBLE  RAY  ANGLES 
FROM  THE  SOURCE.   (THIS 
PROGRAM  USES  ANGLES  FROM 
3  DEGREES  FROM  THE  VERT- 
ICAL TO  177  DEGREES  FROM 
THE  VERTICAL,  IN  3  DEGREE 
INCREMENTS ) 

USE  LINEAR  INTERPOLATION  TO 
CALCULATE  THE  f ( DIFFERENCE ) 
FOR  EACH  INCREMENT  OF  RANGE 
FOR  A  GIVEN  DEPTH 

STORE  ANSWERS  IN  RANDOM 
ACCESS  FILE 

CONTINUE  TO  LOOP  THROUGH 
DEPTH  INCREMENTS  UNTIL 
DONE 


ON 


ERROR 


-Ks) 


RETURN 
TO  MENU 
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7.   Display  Answer  File  On  CRT  Or  List  On  Printer 


DISPLAY  OUTPUT  MENU  TO 
USER  ALLOWING  FOR  SELECTION 
OF: 

1 .  CRT  SCREEN 

2.  PRINTER 

RETRIEVE  HORIZONTAL  RANGE  VS. 

DEPTH  ANSWER  MATRIX  FROM 

RANDOM  ACCESS  FILE  AND 

OUTPUT  TO  USERS  CHOICE 

FROM  ABOVE 


|  RETURN 
(TO  MENU 


<*> 
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8.   Exit  Program  And  Return  To  Operating  System 


0 


WARN  USER  THAT  THIS 

OPERATION  WILL  DESTROY 

ALL  TEMPORARY  FILES 

SUCH  AS  TEMPERATURE 

AND  SOUND  VELOCITY 

PROFILES 


ASK  USER  IF 
EXIT  IS  STILL 
DESIRED,  YES  CR  NO 


YES 


NO 


RETURN 
HTO  MENU) 


EXIT  TO  OPERATING 
SYSTEM 
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VI.   DEFINITION  of  COMPUTER  PROGRAM  VARIABLES 


ADD  Equation  used  to  compute  raypath  length  in 
meters  for  both  the  direct  and  reflected  paths. 

AFLAG  Branching  flag  used  in  solving  for  answer  matrix 
elements 

ANS         Answer  obtained  from  numeric  integration 

AOUT  Answer  matrix  element  output  to  CRT  screen  or 
printer 

BFLAG  Branching  flag  used  in  solving  for  answer  matrix 
elements 

BNDLOW  Lower  limit  for  evaluating  numeric  integration 
II 

BOUNDHI     Upper  limit  for  evaluating  numeric  integration 

BOUNDLO     Direct  path  lower  limit 

CLAT  Latitude  variation  used  in  sound  velocity 
equation 

DDEEP  First  half  of  depth  effect  used  in  sound 
velocity  equation 

DEEP  Second  half  of  depth  effect  used  in  sound 
velocity  equation 

DEG.  Average  whole  degree  latitude  of  operating  area 
used  to  calculate  CLAT 

DELTA  Line  integral  uniform  divisions  used  in  Romberg 
numeric  integration  routine 

DHI  Direct  path  horizontal  range  upper  bound  used  in 
interpolating  direct  path  length  for  desired 
range  increment 

DIRFRAC     Direct  path  interpolation  fraction 

DIRRAY      Direct  path  interpolated  path  length 
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DLONG 

DLOW 

DPH 
DSHORT 

FANSWER.FL1 
FINAL  (  ) 

GRAD 

HYDEEP 

HYDPTH 

HYSPD 

I 

IFLAG 

II 

IN 

INCR 

INN 

J 

JMP 

K 

M 

METER 


Dirct  path  length  upper  bound  used  in 
interpolating  path  length  for  desired  range 
increment 

Direct  path  hoirzontal  range  lower  bound  used  in 
interpolating  direct  path  length  for  desired 
range  increment 

Depth   increment   used   for   outputing   answer   matrix 

Direct  path  length  lower  bound  used  in 
interpolating  path  length  for  desired  range 
increment 

Random  access  file  used  to  store  the  horizontal 
range  vs.   depth   answer  matrix 

Answer  array  used  to  temporarily  hold  the 
interpolated  values  before  storing  them  in 
random  access   file 

Average  temperature  gradilent  input  by  user 

Hydrophone/receiver  depth  in  feet 

Hydrophone/receiver  depth  in  meters 

Sound  velocity  at  hydrophone/receiver  depth 

Loop  increment 

Loop  increment  used  to  branch  through  one  of 
four  integrated  equations 

Loop  increment 

Loop  increment 

Loop  increment 

Loop  increment 

Counter 

Array  element  pointer 

Counter 

Counter 

Array  element  pointer 
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MSTEP 

MSV 

MTRHI 

MTRLO 

MTRS 

MTRSTEP 

MULT 

N 

NANGL 

NMTR 

NN 

NSTEP 

NT 

NTRA 

NUM 

NUMB 

NUPLMT 

NX   (   ) 

PART 

PIECES 

RAD 
RAYDIR   ( 

RAYREF 


Loop  increment 

Array  element  pointer 

Array  element  pointer 

Array  element  pointer 

Array  element  pointer 

Array  element  pointer 

Random  access  file  pointer 

Array  element  pointer 

Ray  angle  from  source 

Horizontal  range  increment 

Counter 

Ray  angle  from  source 

Array  element  pointer 

Array  element  pointer 

Main  menu  index  selection 

Output  menu  index  selection 

Incremented  source  depth 

Array    of    counters     used    during     integration 
accuracy   testing   loop 

Partitioning      factor      used      in     numerical 
integration  routine 

Incremented    divisor    used    to    partition    depth 
integration  limits 

Conversion  of   ray  angle   from  degrees   to   radians 

)    Direct   ray    path    length    calculated    from    equation 
#12    page   25 

Reflected    path     length    calculated    from     equation 
#13   page   25 
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REFFRAC 
REFLO 

REFRAY 
RHI 

RLONG 

RLOW 

RMTR 
RNGDIR   (  ) 

RNGREF  (  ) 

RSHORT 

SALIN 
SNDHI 
SNDLO 
SNDSPS   (  ) 

SNELL 

SOUND 
SST 
START 
SUM 


Reflected  path  interpolation  fraction 

Reflected  path  lower  limit.  Negative  value  of 
boundlo 

Reflected  path   interpolated  path   length 

Reflected  path  horizontal  range  upper  bound  used 
in  interpolating  reflected  path  length  for 
desired  range 

Reflected  path  length  upper  bound  used  in 
interpolating  path  length  for  desired  range 
increment 

Reflected  path  horizontal  range  lower  bound  used 
in  interpolating  reflected  path  length  for 
desired  range 

Depth   increment    for    interpolation   routine 

Direct  path  horizontal  range  array  calculated 
from    equation    §1  page   21 

Reflected  path  horizontal  range  array  calculated 
from    equation   #9  page   22 

Reflected  path  length  lower  bound  used  in 
interpolating   path   length 

Salinity  value  for  operating  area  input  by  user 

Sound  velocity  at  numeric  integral  upper  bound 

Sound  velocity  at  numeric  integral  lower  bound 

Sound  velocity  profile  array  input  directly  by 
user  or  calculated  from  temperature,  salinity, 
and  latitude  inputs 

Snell  constant  equation  result  calculated  from 
initial  ray  angle  and  depth  of  source 

Sound  velocity  at  incremented  source  depth 

Sea  surface  temperature 

Initial  estimate  of  integral  answer 

Equation  used  to  compute  horizontal  range  in 
meters  for  both  the  direct  and  reflected  paths 
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SVP         Input  sound  velocity  profile  by  user 

TEMP        Degrees  celsius  temperature  divided  by  1 0 

TEST  (  )     Successive  numeric  integration  answer  array  used 
to  check  accuracy  of  answer 

TGTDPTH     Incremented  source  depth 

TOL  Accuracy    desired    for    numerical     integration 

routine  answer 

TOT  Summation    of    partitioned  numerical    ilntegration 

values 

TTEMP  Input     temperature    used    in    calculating     sound 

velocity  profile 

WALK        Integration  routine  smoothing  factor 

X  Incremented  depth  value  for  each  partition 
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APPENDIX  A 

Urick  [Ref.  1:  p.  124]  briefly  discusses  the  intensity  of  sound 
reaching  a  receiver  via  a  direct  path  and  a  surface  reflected  path.  The 
following  derives  the  ratio  of  the  total  intensity  of  the  signal  received 
at  a  point  when  multipath  Lloyd's  mirror  interference  is  present,  to  the 
total  intensity  resulting  only  from  the  direct  path. 


SURFACE 


IT- 


DEPTH 

METERS 


it 


MUST  SHOW: 


REFLECTED 


'DIRECT 


RANGE 
METERS 


S0URCE       <I  >  +  <I  > 

Intensity  Ratio  =  /T  r  = 77-r 

^0        0 


Rr 


%  ■ '  *  ■  (s!  -  >%  -*■  > 

WHERE:  <I>,  the  time  averaged  intensity  and  I,  the  instantaneous 
intensity  are  defined  as  follows: 


(1)  <I>  (l/T) 


1 

/' 


(2) 


T  =  1  period 


(P0(r)  +  PR(r))2 
I  *  — or^ J  Pn(D 


-  ,  ,  QVI ,  and  P^ir)   are  both  functions  of  time 
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GIVEN: 


M    Deceived  ::  pdirect  +  preflected 


po  yP0 

=    n~    COS    U)t    -   5 COS    [u>(t-T)] 

KD  KR 

where:     y  =  surface  pressure  reflection  coefficient 

t  =  time  delay  caused  by  the  difference   in  the  two  path 
lengths 
2 
2 


2   P 
u   pc   2pcR 


D 


where:  P  =  RMS  pressure  amplitude  of  the  direct  path  ray 

NOTE:  It  is  assumed  that  the  transmission  loss  is  due  only  to  spherical 

spreading  and  that  absorption  in  the  medium  is  negligible. 
SOLUTION: 

P0  -mPo 

PDIR  =  R"~  C0S  ^  PREF  =  R COs[u>(t-r)] 

U  K 


PRECEIVED  =  P0IRECT  +  PREFLECTED 


PRECEIVED  =  P0 


TT-  COS(a)t)     -    jj-      COS    (uj(t-T)) 


2 
2        P 

I  =  £_  =  _° 

pC    ~2pC 


2  2 

COS    (gjt)  u       „nrZ(     I*.       \\ 

^ L   +  — ~    COS    (oj(t-x)) 


.2a. 


RRRD 


COS    (ut)    COS    Co>(t-x)] 
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■+/ 


SINCE:     <I>  =  4     /I     dt 


AND       a)  =  2irf  =  y- 


2tt/oj 

I     dt 


<I>  = 


u  p 


2irpC 


2ir/o)  «       2 

-^      /      cos2(wt)  dt  +  Uj    [ 


2tt/oj 

\    /        cos2(o3(t-x))  dt 


2u 
RDRR 


2tt/, 


/  COS(o)t 


)  cos(u)(t-T))  dt 


Solving  the  three  integrals  separately 


2tt/w 
1        /*  2 

—      /        COS    (cot 

RD   ± 


)   dt  =  ^  [|+  ^  sin(2  ut)] 
RD 


=  Kt-+T"  sin(4Tr)] 
n2     (i)       4co  v      ' 


2tt/w 
0 
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2tt/ 


U) 


%2  cosc(o)(t-T);    dt 

R 


2 

ii_ 

R2 
RR 


|  +  ^  sin(2o,(t-x)) 


2w/ 


CO 


_  y 


i  +  ^  (sin(4*  -  2UT)-s1n(-2UT)) 


RRRD 


2tt/o 


cos(wt)  cos(w(t-x))  dt 


_  2y 


RRRD 


/" 


CO 


cos  (wt)  cos(cox)+sin(cot)   cos(u>t)  sin(a)x)  dt 


_  2u  cos(cj 


RRRD 


2ir/w  2tt/co 

^      /"         COS2(cot)    dt  +   2HrS™(u>t)       /*         s1n(ajt)    cos(ut) 


dt 


_  2y  cos(gjx )     ft   ,    1     e4„/o  *\ 
3       rr    —      (t  +  XT  sin(2wt) 

KRKD 


■2       4oj 


2ir/a: 


0 


+  u?lfi(sin2u)) 


VR,XD 


27r/a, 
0 


_  2u  cos  jmt  )    /jr»  +  2y  sin((jax )    /qx 
RRRD  a)  RRRD 


(^V~  cos(wx  ) 
KRRD 
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Assembling  the  answers  to  the  three  integrals  yields: 
2 


U)P 


<I>  = 


0 


2lTpC 

n2     r 


tt  \x   tt       2 iry        „_. 

—  '^2    "  Z^  cos^> 


wRD       wRR 


<I>  = 


<I>  = 


r0 
2pc 

P2 


2pcRp 
NOTE: 


5"  +    p  "  p  d     cos(WT ) 
Rjj       R^       KRKD 


M  RD       2yRD 
1   +  — s p—  cos(WT) 

RR  R 


<I0>  = 


2PcR 


<I>  =  <IQ> 


Finally 


u  RD       2pRD 


1  +  — 2 R —  COS(ujt) 


R 


This  result  shows  the  dependence  of  the  intensity  at  the  receiver  on 
the  difference  in  path  lengths  RnjocrT  and  rrfflectfd  anc'  tne  t1me  delay, 
t,   caused  by  this  difference.     This   intensity  ratio  is  a  minimum  (corre- 
sponding to  a  null   on  the  Lloyd's  mirror  interference  pattern)  when 
cos  ojt  =  1  or  got  =  integer  multiples  of  2tt. 
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COMPUTER  OUTPUT 

A  sample  of  the  actual  computer  output  is  listed  below. 
Data  values  are  frequency  measurements  Af  (Hertz)  between 
successive  nulls  on  the  Lloyd's  mirror  pattern  (see  Figure 
II-l). 

Given  a  particular  source  depth  and  a  Af,  the  source  range 

can  be  determined  from  the  corresponding  output  column.   The 

source  depth  can  be  determined  from  the  appropriate  output  row 

if  the  horizontal  range  to  the  source  from  the  receiver  and 

the  Af  are  known. 

Group  of  Ranges  From  10  to  70  Meters 

Depth  of  10  Meters 

80.6     94.8    114.1    136.0    159.8    184.4    210.1 

Depth  of  20  Meters 

41.9     50.2     60.1     71.0     82.7     94.8    107.3 

Depth  of  30  Meters 

31.6     37.3     43.8     50.8     58.3     66.1     74.3 

Depth  of  40  Meters 

29.2     32.5     36.8     41.7     46.9     52.5     58.3 

Depth  of  50  Meters 

28.5     30.6     33.5     36.9     40.8     44.9     49.4 
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COMPUTER    PROGRAM 

10     DEFINT     I,J,K,L,M,N 

20    DIM    NX(16)  ,TEST(1  36 ) , SNDSPD(  201 ) ,RNGDIR( 1 80) 

30   DIM   RNGREF(180),RAYDIR(180),RAYREF(1  80),FINAL(1  01  ) 

40    TOL    =    .05 

50    OPEN    "R",#1  ,  "FANSWER.FL1",4 

60  FIELD   #1  ,    4  AS    F1  $ 

100    PRINT    CHR$(26):     PRINT    "SELECT    OPERATION":  PRINT 

110   PRINT   1  /'INITIALIZE  ANSWER  FILE, DO  AT  THE    BEGINNING  OF 

EACH   NEW    BT    INPUT" 

120    PRINT    2, "INPUT    BT    INFO    IN    5    METER    INCREMENTS" 

130   PRINT   3 /'INPUT  SEA  SURFACE    TEMP   AND   GRADIENT  FOR    FIRST 

200    METERS" 

135     PRINT     4, "INPUT     SOUND    VELOCITY    PROFILE     IN    5     METER 

INCREMENTS" 

140     PRINT     5/'CALCULATE     LLOYD'S    MIRROR     HORIZONTAL    RANGE     .VS. 

DEPTH    CHART" 

150    PRINT    6  /'OUTPUT    CHART   TO    PRINTER" 

160    PRINT    7, "EXIT    PROGRAM" 

170    PRINT:PRINT:INPUT    "ENTER    THE    NUMBER    OF    THE    OPERATION 

DESIRED         ",NUM 

180    IF    (NUM<1)    OR    (NUM>7)    THEN  PRINT   "IMPROPER   INPUT":GOTO 

100 

190    ON    NUM    GOSUB    1000,2000,5000,5500,6000,9600,10000 

200    GOTO    100 

1000   REM      INITIALIZE  ANSWER    FILE 

1010    PRINT    CHR$(  26): PRINT    "ARE    YOU    SURE    YOU    WANT    TO    ZERO    OUT 

THE    ANSWER    FILE" 

10  30    INPUT   "ENTER  Y  FOR   YES,    N  FOR  NO         ",A$ 

1040    IF    A$<>"Y"    THEN    RETURN 

1050    CLOSE    #1 

1060    KILL    "FANSWER.FL1" 

10  70    OPEN    "R",#1  ,"FANSWER.FL1",4 

1080  FIELD   #1  ,    4  AS  F1  $ 

1090    RETURN 

2000    REM    INPUT    BT    INFO    IN    5    METER    INCREMENTS,  THEN    CALCULATE 

SOUND   VELOCITY 

2010    REM  USING   ALAN   R.    COPPENS    EQUATION    PRESENTED    IN    J. 

ACOUS.    SOC.    OF   AMER. 

2015  REM  VOLUME     69,     PAGE     862,      1981     AND     PUBLISHED     IN 
FUNDAMENTALS   OF    ACOUSTICS 

2016  REM  THIRD     EDITION,      BY     KINSLER    FREY     COPPENS     AND 
SANDERS 

2020    GOSUB    3010:     GOSUB    4000:PRINT    CHR$(26) 

2030    PRINT    "ENTER    BT    INFO    IN    5    METER    INCREMENTS    STARTING 

WITH    THE    SURFACE" 

2040    FOR    I    =    0    TO    200    STEP    5 

2050     PRINT:PRINT     I     "M  ETERS" :  INPUT    "TEMP     IN     CELSIUS     = 

",TTEMP:PRINT 

2060  PRINT  "YOU  INPUT  "TTEMP"  DEGREES  CELSIUS": 
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INPUT    "IS     THIS     CORRECT     ?     ENTER    Y    FOR    YES       N     FOR    NO 
";B$ 

2070    IF   B$o"Y"    THEN   PRINT    "REENTER   TEMP":GOTO    2050 
2075   TEMP   =  TTEMP   /    10! 
2080  METER  =   I 

2100    DEEP     =    (16.23    +.253*TEMP)*    CLAT    *    CSNG(  I)  /1 000!      +  (.213- 
.1*TEMP)     * 

(CLAT    *    (CSNG(I)/1000!))22 
2105     DDEEP     =      (.016      +     .0002     *      (SALIN-35))      *      (SALIN- 
3  5 ) *TEMP*CLAT*CSNG ( I ) / 1  0  0  0 ! 

2110  SNDSPD(METER)  =1449.05    +    45.7    *    TEMP    -    5.21     *    (TEMP22)    + 
.23   *    (TEMP23)    + 

(1  . 333-.  12 6*TEMP+. 00 9*(TEMP22 ) ) *( SALIN- 
35  )+(DEEP+DDEEP) 

21 11  IF    I    =    0    GOTO    2120 

2112  SLOW    =   SNDSPD(    I   -    5   ) 
21  1  3    SHI    =    SNDSPD(    I    ) 

21 1 4  FOR   J    =    1    TO    4 

2115  SNDSPD(  I  -  5  +  J  )  =  SLOW  +  (  J  /  5  )  *  (  SHI  -  SLOW  ) 

2116  NEXT  J 
2120  NEXT  I 
21  30  RETURN 

3000  REM   ENTRY  OF  OPAREA  LATITUDE 

3010  PRINT  CHR$(  26): PRINT  "ENTER  AVERAGE  DEGREES  LATITUDE  OF 

OPAREA" 

3020  INPUT  "IN  WHOLE  DEGREES     ",DEG 

3030  PRINT:PRINT  "YOU  ENTERED  "DEG"  DEGREES,  IS  THIS  CORRECT 

?" 

3040  INPUT  "ENTER  Y  FOR  YES   N  FOR  NO    ",CC$ 

3050  IF  CC$o"Y"  THEN  PRINT  "REENTER  LATITUDE  ":GOTO  3010 

3060  CLAT  =   1  -  .0026  *  (COS(  DEG  *  3.14159  /  90!)) 

3070  RETURN 

4000      REM         ENTRY     OF      SALINITY     VALUE     FOR      SOUND     SPEED 

COMPUTATIONS 

4010    PRINT    CHR$(26) 

40  20    PRINT   "ENTER  THE    SALINITY    VALUE  AS   PARTS   PER  THOUSAND 

FOR    THE    OPAREA" 

4030     INPUT    "IF    UNKNOWN,     ENTER    A    VALUE    OF    35,     SALINITY    = 

",SALIN 

4040    PRINT:PRINT    "THE    SALINITY    LEVEL    IS    "SALIN"    PARTS    PER 

THOUSAND" 

40  50     INPUT    "IS     THIS     CORRECT    ?     ENTER    Y    FOR    YES       N    FOR    NO 

",BB$ 

4060    IF  BB$    <>   "Y"   GOTO   4010 

4070    RETURN 

5000    REM    ENTER    SEA    SURFACE    TEMP    AND   AVERAGE    GRADIENT    FOR 

FIRST   200    METERS 

5010    PRINT    CHR$(26) 

50  20    INPUT    "ENTER    THE    SEA    SURFACE    TEMPERATURE    IN    DEGREES 

CELSIUS  ",SST 

5030  PRINT  .-PRINT  "YOU  ENTERED  "SST"  DEGREES  CELSIUS  FOR  THE 
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SEA    SURFACE    TEMP" 

50  40    INPUT    "IS    THIS    CORRECT    ?       ENTER    Y    FOR    YES       N    FOR    NO 

",C$ 

5050  IF  C$  <>  "Y"  THEN  PRINT  "REENTER  SEA  SURFACE  TEMP": GOTO 

5020 

5060  PRINT 

5070      PRINT      "ENTER     THE      TEMPERATURE     GRADIENT      (  DE  G 

CELSIUS /METER)" 

5075  INPUT  "FOR  THE  FIRST  200  METERS     ",GRAD: PRINT 

5080  PRINT  "YOU  ENTERED  A  TEMPERATURE  GRADIENT  OF  "GRAD" 

DEGREES  CELSIUS /METER" 

50  90  INPUT  "IS  THIS  CORRECT  ?   ENTER  Y  FOR  YES   N  FOR  NO 
",D$ 

5100  IF  D$  <>  "Y"  THEN  PRINT  "REENTER  GRADIENT": GOTO  5070 

5110  GOSUB  3010 

5120  GOSUB  4000 

5125  PRINT  CHR$(26)  .-PRINT  "PLEASE  WAIT,  CALCULATING  SOUND 

VELOCITY  VALUES" 

5130  FOR  I  =  0  TO  200 

5140  TEMP  =  (SST  +  (CSNG(I)  *  GRAD))  /  10! 

5150  DEEP  =  (16.23  +  .253*TEMP)*  CLAT  *  CSNG(I)  / 1 0  00 1   + 

( .213-.1*TEMP) 

*  (CLAT  *  (CSNG(I)  /1000!  ))22 

5155    DDEEP    =     (.016    +    .0002    *     (SALIN    -    35))    *     (SALIN-35)     * 

TEMP   *  CLAT  * 

CSNG(I)    /    10  00! 
5170       SNDSPD(I)    =    1449.05    +    45.7    *  TEMP  -    5.21    *    (TEMP22)    + 
.23   *    (TEMP23    ) 

+(1.333       -     .126      *      TEMP      +      .00 9*TE MP a 2 ) * ( SALIN      - 
35)+DEEP+DDEEP 
5180    NEXT    I 

51  90    RETURN 

5510  REM  ENTRY  OF  SOUND  VELOCITY  PROFILE 

5525  PRINT  CHR$(26) 

5530     PRINT     "ENTER     SOUND     VELOCITY     PROFILE     IN     5     METER 

INCREMENTS    " 

55  40     PRINT    "STARTING    WITH    THE    SURFACE     SOUND    VELOCITY    IN 

METERS /SECOND" 

5550    FOR    II    =    0    TO    200    STEP    5 

5560       PRINT:PRINT       II       "METERS"       :INPUT       "       SOUND 

VELOCITY(METERS/SEC)=       ",SVP 

5570    PRINT: PRINT    "YOU    INPUT    "SVP"    METERS/SEC" 

55  80     INPUT    "IS     THIS     CORRECT    ?     ENTER    Y    FOR     YES     N    FOR    NO 

";SS$ 

55  90    IF    SS$<>"Y"    THEN    PRINT   "REENTER    SOUND   VELOCITY":GOTO 

5560 

5610  SNDSPD(II)    =    SVP 

5611  IF    II    =    0    GOTO    5620 

5612  SLOW    =    SNDSPD(    II    -    5    ) 

5613  SHI   =  SNDSPD(    II    ) 
561  4    FOR   JJ    =    1    TO    4 
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5615  SNDSPD(     II    -    5    +    JJ    )     =    SLOW    +     (    JJ    /    5    )     *    (     SHI    - 
SLOW    ) 

5616  NEXT  JJ 
5620  NEXT  II 
5630    RETURN 

60  00     REM    CALCULATE     LLOYD'S     MIRROR     HORIZONTAL    RANGE     .VS. 
DEPTH    CHART 

6005  PRINT:PRINT:PRINT    "WILL    HYDROPHONE   DEPTH   BE    INPUT   IN 
FEET    OR   METERS    ?" 

6006  INPUT   "ENTER    F    FOR    FEET    OR   M    FOR    METERS  ",     AHYD$ 

6007  IF  AHYD$   =  "F"   GOTO  6010   ELSE  IF  AHYD$   =  "M"   GOTO  6052 

6008  PRINT: PRINT: PRINT    "INPUT    ERROR"    :GOTO    60  05 

6010    PRINT    CHR $( 26 ): INPUT    "ENTER    HYDROPHONE    DEPTH    IN    FEET 

"fHYDEEP 

6020  PRINT :PRINT  "YOU  ENTERED  A  HYDROPHONE  DEPTH  OF  "HYDEEP" 

FEET" 

60  30     INPUT    "IS     THIS     CORRECT    ?     ENTER    Y    FOR    YES       N    FOR    NO 

",AA$ 

6040    IF  AA$    <>   "Y"   GOTO  6010 

6050  HYDPTH    =    (HYDEEP    *    .3048) 

6051  GOTO    6060 

6052  PRINT    CHR$(  26):  INPUT    "ENTER    HYDROPHONE    DEPTH    IN    METERS 
", HYDPTH 

6053  PRINT:PRINT     "YOU     ENTERED     A     HYDROPHONE    DEPTH     OF 
"HYDPTH"     METERS" 

60  54    INPUT    "    IS    THIS    CORRECT       ?       ENTER    Y    FOR    YES    N    FOR   NO 
",BHY$ 

6055   IF  BHY$    <>   "Y"  GOTO  6052 
6060    HYSPD    =    SNDSPD(    HYDPTH    ) 

61  70   FOR  NUPLMT  =   1  0   TO  200  STEP  1  0 
61 75    PRINT   "NUPLMT="    NUPLMT 

6180    TGTDPTH    =    CSNG(    NUPLMT    ) 

6185    IF    HYDPTH    =    TGTDPTH    THEN    HYDPTH    =    HYDPTH    -    3 

6190    IF   TGTDPTH    >    HYDPTH   GOTO    6230 

6200    BOUNDLO    =    TGTDPTH 

6210    BOUNDHI    =    HYDPTH 

6220    GOTO    6250 

6230    BOUNDLO    =    HYDPTH 

6240    BOUNDHI    =    TGTDPTH 

6250    REFLO    =    -(BOUNDLO) 

6360    FOR   NSTEP    =    3    TO    1 77    STEP    3 

6365    PRINT    "NSTEP="    NSTEP 

6370  RAD  =  CSNG(NSTEP)   *  3.14159   /    180! 

6380    MTRSTEP   =    TGTDPTH 

6390    SOUND    =    SNDSPD(MTRSTEP) 

6400    MTRHI   =    BOUNDHI 

6410    SNDHI    =    SNDSPD(MTRHI) 

6420    MTRLO  =    BOUNDLO 

6430    SNDLO    =    SNDSPD(MTRLO) 

64  40    SNELL    =   SIN(    RAD   )    /    SOUND 

6480    FOR    IFLAG   =    1    TO    4 
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7000   REM   ROMBERG   NUMERICAL   INTEGRATION    SUBROUTINE 
7015   START  =  01 
7020    PIECES    =    1! 
7030    NX(1 )    =    1 

7040  IF     (IFLAG  =  2)      OR      (IFLAG  =  4)      THEN     BNDLOW=REFLO     ELSE 
BNDLOW=BOUNDLO 

7041  IF    (IFLAG    =   2)    OR   (IFLAG  =   4)    GOTO  7045 

7042  IF    (TGTDPTH    >    HYDPTH)    AND    (NSTEP    >     90)    GOTO    8105 

7043  IF    (TGTDPTH    <    HYDPTH)    AND    (NSTEP    <     90)    GOTO    8105 

7044  GOTO    7050 

7045  IF   NSTEP    >    90    GOTO    8105 

7050  DELTA    =    ( BOUNDHI    -    BNDLOW)    /    PIECES 

7051  DENOM1     =    (SNELL) 22    *    (SNDL0)22 

7052  DENOM2    =    (SNELL) 22    *    (SNDHI)22 

7053  IF    DENOM1     =    1!    THEN    DENOM1     =   .99  9   :IF  DENOM2    =    1!    THEN 
DENOM2    =    .999 

7054  IF    (    DENOM1    >    1!    )    OR    (    DENOM2    >    1!    )    GOTO  8105 
7060    IF    (IFLAG=1)    OR    (IFLAG=2)     GOTO    7085 

7065    FOR    IN    =    1    TO    2 

70  68    IF   IN   =  1    THEN  DENOM   =   DENOM1    ELSE   DENOM  =  DENOM2 

7070    START    =    START    +     ((1!)/    SQR(1  !     -    DENOM)) 

7075    NEXT    IN 

7080  GOTO  7100 

7085  FOR  INN  =  1  TO  2 

7088  IF   INN  =   1    THEN  DENOM   =  DENOM1    ELSE  DENOM  =  DENOM2 

7089  IF    INN    =    1     THEN    SPEED    =    SNDLO    ELSE    SPEED    =    SNDHI 

7090  START    =    START    +    ((SNELL    *    SPEED)     /    SQR(1!     -    DENOM)) 

70  95    NEXT    INN 

71  00    TOT  =  START  /   2! 

71  1  0   TEST(  1  )    =   DELTA  *   START    /    2  ! 

7120    N    =    1 

7130    NN    =    2 

7140    N    =    N    +    1 

7150    WALK    =    41 

71 60    NX(N)    =  NN 

7170    MSTEP    =    CINT(    PIECES    ) 

7180    PIECES    =   PIECES    *    2! 

7190    DELTA    =    (BOUNDHI    -    BNDLOW)    /    PIECES 

7500   REM  COMPUTE   TRAPEZOIDAL   SUM 

7510  FOR  INCR  =   1    TO  MSTEP 

7520   PART  =   CSNG(INCR)   *   2i    -    1! 

7530    X    =    BNDLOW    +    (DELTA    *    PART) 

7540    MTRS    =   CINT    (    ABS    (X)    ) 

7545   IF  MTRS    >  200   THEN   MTRS   =   200 

7550    SPEED    =    SNDSPD(    MTRS    ) 

7552    DENOM3    =    (SNELL) 22    *    ( SPEED) 22 

7554  IF   DENOM3    =    11     THEN   DENOM3    =    .999 

7555  IF    DENOM3    >    1!    GOTO    8105 

7570    IF    (IFLAG=1)    OR    (IFLAG=2)    GOTO    7600 
7580   TOT  =    TOT  +    (    (1!)    /    SQR(1!    -    DENOM3)) 
7590    GOTO    7610 
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7600   TOT   =  TOT  +    (    (SNELL   *   SPEED)    /    SQR(1l    -    DENOM3)) 

7610    NEXT    INCR 

7620    TEST(    NN    )    =   TOT   *  DELTA 

7630    K    =   N    -    1 

7640    NTRA  =    NX(    K    ) 

76  50   FOR   M   =   1    TO   K 

7660    J  =   NN    +  M 

7670   NT   =  NTRA  +    M    -    1 

7680    TEST(J)     =    (     TEST(J-1)     *    WALK    -    TEST(NT))     /    (WALK-1!) 

76  90  WALK  =  WALK  *   4! 

7700  NEXT  M 

7705    IF    (IFLAG    =    1    OR    IFLAG    =   2)    AND    (TEST(J)    >    3000!)    GOTO 

8105 

7710    IF   N    <    5    GOTO    7760 

7720    IF   TEST(    NN  +    1     )    =   0!    GOTO   7750 

7730     IF    ABS(TEST(NTRA+1 )-TEST(NN     +     1))     <=     (ABS(TEST(NN     + 

1)*TOL))     GOTO    8000 

7740    IF   ABS(    TEST(NN-1  )    -    TEST(J))    <=    (ABS(    TEST(J)    *    TOL)) 

GOTO    8000 

7750    IF   N    >    15    GOTO    9500 

7760    NN    =   J    +    1 

7770  GOTO  7140 

8000  ANS  =  TEST(J) 

8030    ON    IFLAG   GOTO    8040,8060,8080,8100 

8040    RNGDIR(NSTEP)    =    ANS 

8050    GOTO    8105 

8060    RNGREF(NSTEP)    =    ANS 

8070    GOTO    8105 

8080    RAYDIR(NSTEP)    =    ANS 

8090    GOTO    8105 

8100    RAYREF(NSTEP)    =    ANS 

8105    NEXT    IFLAG 

81  10    NEXT   NSTEP 

8115    REM       USE    LINEAR    INTERPOLATION  TO   CALCULATE  ANSWER    CHART 

VALUES    IN   RANGE 

8117    REM  AND    DEPTH   INCREMENTS    OF    10    METERS 

8120   FOR  NMTR  =1    TO   100 

8130    RMTR  =  CSNG(  NMTR  )    *   10! 

8135   DIRRAY   =   0!    :    REFRAY    =   0! 

81  40  AFLAG  =   0! 

81 50   BFLAG  =   0! 

8160  IF   TGTDPTH    <    HYDPTH   GOTO   8179 

8161  FOR   NANGL    =    3    TO    87    STEP    3 

8162  IF    RNGDIR(    NANGL    )     <=    RMTR   GOTO    8171 

8163  DLOW    =    RNGDIR(    NANGL    -    3     ):DHI    =    RNGDIR(    NANGL    ) 

8164  IF   DLOW    =    DHI    THEN   DIRFRAC    =    0    :GOTO    8166 

8165  DIRFRAC    =    (RMTR    -    DLOW)    /     (DHI    -    DLOW) 
81  66   IF  DLOW   >    0  GOTO  8168 

8167  DSHORT    =    BOUNDHI    -    BOUNDLO    :GOTO    8169 

8168  DSHORT    =    RAYDIR(    NANGL    -    3    ) 

8169  DLONG    =    RAYDIR(    NANGL    ) 
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81  70 

8210 
81  71 
ELSE 
8172 
8173 
8174 
81  75 


DIRRAY    =    DSHORT    +    (    DLONG    -    DSHORT    )    *   DIRFRAC    :     GOTO 
=    87     )    AND    (    RNGDIR(    87    )     >    0    )    GOTO    8172 


IF    (     NANGL 
GOTO    8177 

DLOW    =   RNGDIR(    87    ) 

IF    DLOW    =    RMTR    THEN 

DIRFRAC      =    (    RMTR    - 

DSHORT     =     RAYDIR( 

(BOUNDHI    -    BOUNDLO)22     ) 

8176    DIRRAY    =    DSHORT    + 

8210 

8177 

8178 

81  79 

8180 

81  81 

ELSE 

8182 

8183 

81  84 


DIRRAY 
DLOW    ) 
NANGL 


RAYDIR( 
(    1000    - 
:      CLONG 


87    ) 

DLOW    ) 
=     SQR( 


GOTO    821 0 


100022     + 


(    CLONG    -    DSHORT    )    *   DIRFRAC    :    GOTO 


NEXT  NANGL 

GOTO  8410 

FOR   MANGL  =   93   TO  177   STEP 

IF   RMTR    >=    RNGDIR(    MANGL    ) 


3 
GOTO 


8181    ELSE    GOTO    81 95 


=    93     )    AND    (    RNGDIR(    93     )     >    0    )    GOTO    8182 


IF    (    MANGL 
GOTO    8186 

DLOW    =   RNGDIR(    93    ) 
IF    DLOW    =   RMTR   THEN    DIRRAY 
DIRFRAC  =    (    RMTR   -    DLOW    ) 
RAYDIR(     93     )    : 

BLONG  =  SQR    (   100022   +   (    BOUNDHI   -    BOUNDLO   )22   ) 


=    RAYDIR( 
/    (    1000    - 


93    )    : 
DLOW 


GOTO    8210 
)    :    DSHORT    = 


81  85 

8210 

8186 

8187 

GOTO 

8188 

81  89 

81  90 

8191 

8192 

81  94 

8210 

81  95 

81  96 

8197 

81  98 

8210 

8199 

8200 

8205 

8210 

8211 

8212 

8213 

8214 

8215 

8216 

8217 

8219 


DIRRAY=    DSHORT    +    (    BLONG    -    DSHORT    )    *    DIRFRAC    :     GOTO 


=    RNGDIR(    MANGL    )    : DHI    =    RNGDIR(    MANGL    -    3    ) 


DLOW 

IF    ( 

8195 

IF   DLOW    = 

DIRFRAC  = 

IF  DLOW    > 


DLOW    <=  RMTR    )    AND    (    DHI    >=  RMTR    )    GOTO    8188    ELSE 


DHI    THEN   DIRFRAC    =    0    :    GOTO    81 90 
(  RMTR   -    DLOW    )    /    (   DHI   -    DLOW    ) 
0  GOTO  81  92 
DSHORT    =    BOUNDHI    -    BOUNDLO   :     GOTO    8193 
DSHORT    =    RAYDIR(     MANGL    ):    DLONG    =    RAYDIR(    MANGL    - 
DIRRAY    =   DSHORT   +    (    DLONG    -    DSHORT    )    *    DIRFRAC    : 


IF    (    MANGL   <>    177   )    GOTO  8199 
DLOW    =   0   :    DHI   =  RNGDIR(    177    ) 
DSHORT    =    BOUNDHI    -    BOUNDLO    : 
DIRRAY    =   DSHORT   +    (    DLONG    - 


DIRFRAC    =  RMTR   / 
DLONG    =    RAYDIR(     177 
DSHORT    )    *    DIRFRAC    : 


3) 

GOTO 


DHI 


GOTO 


NEXT   MANGL 

GOTO    8410 

IF    DHI    =    DLOW    THEN 

FOR   NREF    =    3    TO    87 


=    0 


GOTO    8211 


) 


DIRFRAC 

STEP  3 
IF  RNGREF(  NREF  )  <=  RMTR  GOTO  8220 
RLOW  =  RNGREF(  NREF  -  3  )  :  RHI  =  RNGREF(  NREF 
IF  RLOW  =  RHI  THEN  REFFRAC  =  0  :  GOTO  8215 
REFFRAC  =  (  RMTR  -  RLOW  )  /  (  RHI  -  RLOW  ) 
IF  RLOW   >    0  GOTO  821  7 

RSHORT    =    BOUNDHI    +    BOUNDLO    :     GOTO    8218 
RSHORT    =    RAYREF(     NREF    -    3    ):    RLONG    =    RAYREF(NREF) 
REFRAY    =    RSHORT    +    (    RLONG    -    RSHORT    )    *    REFFRAC 


GOTO 
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8400 

8220  IF    (    NREF    =    87    )    AND    (RNGREF(87)    >    0    )    GOTO    82  21    ELSE 
GOTO    8226 

8221  RLOW    =   RNGREF(    87    ) 

8222  IF   RLOW    =    RMTR   THEN   REFRAY    =    RAYREF(    87    ):     GOTO    8400 

8223  REFFRAC    =    (    RMTR    -   RLOW    )    /    (    1000    -   RLOW    ) 

82  2  4     RSHORT     =    RAYREF(     NREF     )     :     ALONG     =    SQR     (     100022    +     ( 
BOUNDHI    +    BOUNDLO    )22    ) 

82  25    REFRAY    =  RSHORT   +       (    ALONG   -    RSHORT   )    *  REFFRAC   :    GOTO 
8400 

8226  NEXT   NREF 

8227  GOTO    8410 

83  70    REM      OUTPUT    FINAL   ANSWERS  FOR  EACH  DEPTH  AND  RANGE   TO 
DISK    FILE 

8385    LPRINT    "HYSPD="    HYSPD 

8400    FINAL(    NMTR    )    =    HYSPD    /     (    REFRAY    -    DIRRAY    ) 

8410  NEXT   NMTR 

8420    MULT    =    ( TGTDPTH    ±    10)    -    1 

8430  FOR  JMP  =   1    TO   100 

8440    LSET    F1  $    =    MKS$(  FINAL(  JMP) ) 

8450  RECNUM%  =  JMP  +   (MULT  *   1  00) 

8460    PUT    #1,    RECNUM% 

8470  NEXT   JMP 

8480  FOR  NUL  =  0  TO  100 

8490  FINAL(NUL)  =  0» 

8500  NEXT  NUL 

8510  FOR  NUL  =  0  TO  1 35 

8520  TEST(NUL)  =  0! 

8530  NEXT  NUL 

8540  FOR  NUL  =  0  TO  1 5 

8550  NX(1 5)  =  0! 

8560  NEXT  NUL 

8570  FOR  NUL  =  0  TO  1 79 

85  80RNGDIR(NUL)=0!:RNGREF(NUL)=0:RAYDIR(NUL)  =  0: 

RAYREF(NUL)=0 
8620    NEXT   NUL 
9000  NEXT   NUPLMT 
9010    RETURN 

9500    PRINT: PRINT: PRINT    "INPUT    DATA    HAS    CAUSED    AN    ERROR" 
9510    GOTO    100 

9600    REM    PRINT   OUT    GRAPH    OF    HORIZONTAL    RANGE    .vs.     DEPTH 
9610    PRINT    CHR$(26) 

9620    PRINT :PRINT:PRINT    "    SELECT    OUTPUT    DEVICE    ":     PRINT 
9630    PRINT    1,"   VIDIO    SCREEN    ":    PRINT 
9640    PRINT    2,"    PRINTER    ":     PRINT 

9650     INPUT    "    ENTER    THE    NUMBER    FOR    OUTPUT    DEVICE    DESIRED 
",NUMB 

9660    IF    (    NUMB    <    1    OR   NUMB    >    2    )   THEN  PRINT  "IMPROPER   INPUT 
":GOTO    9620 
9665    PRINT: PRINT 
96  70    PRINT    "         THE    OUTPUT    MATRIX    IS    PRESENTED    IN    10    METER 
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INCREMENTS    FOR   BOTH   " 

9680    PRINT   "   HORIZONTAL    RANGES    AND  FOR    DEPTHS  OF    10    METERS 

TO  2  00    METERS.    " 

96  90    PRINT    "         HORIZONTAL    RANGES    FROM    THE   RECEIVER    TO    THE 
TARGET    ,    IN   " 

9700    PRINT    "    GROUPS    OF    100    METERS,     INCREASE    FROM    LEFT    TO 

RIGHT,    " 

9710    PRINT    "   ie.    10    TO    100    METER   RANGES    WILL    BE    DISPLAYED 

FIRST,     FOLLOWED" 

9720     PRINT     "    BY     110     TO    200    METERS    AND    SO    ON    UP    TO    1000 

METERS . " 

97  30    PRINT   "         DEPTHS    OF    0   TO   200    METERS    WILL   BE   DISPLAYED 
FOR  EACH  GROUP   OF  " 

9740     PRINT    "    RANGES,     WITH    THE    FIRST    ROW    REPRESENTING    10 

METERS,     THE   " 

9750  PRINT  "  NEXT  ROW  20  METERS,  ON  DOWN  THE  PAGE  UNTIL  ROW 

20  IS  " 

9760  PRINT  "  PRINTED  FOR  A  DEPTH  OF  200  METERS.":  PRINT:  PRINT 

9770  FOR  I  =  0  TO  9 

9780  THIS  =  10  +  100  *  I 

9790  THAT  =  1 00  +  1 00  *  I 

9800  IF  NUMB  =  1  THEN  PRINT: 

PRINT       "GROUP       OF       RANGES       FROM       "THIS"       TO 
"THAT"METERS":GOTO    9  82  0 
9805     LPRINT: LPRINT: 

9810    LPRINT:LPRINT    "GROUP    OF    RANGES    FROM    "THIS"    TO    "THAT" 
METERS" 

9820  FOR  J  =  0  TO  19 
9830  DPH  =  10  +  10  *  J 

9840      IF      NUMB      =      1      THEN      PRINT:PRINT     "DEPTH     OF     "DPH" 
METERS": GOTO    9860 

9850    LPRINT:LPRINT    "DEPTH    OF    "DPH"    METERS" 
9860    FOR   K    =    1    TO    1 0 

9870    RECNUM%   =    K  +    100   *   J  +   10    *   I 
9880    GET    #1,RECNUM% 

98  90    AOUT    =  CVS(    F1  $    ) 
9900   IF  NUMB  =   2   GOTO   9930 

9910    PRINT    USING    "####.#       ";AOUT; 

9920    GOTO    9940 

9930    LPRINT    USING    "####.#       ";AOUT; 

9940    NEXT   K 

9950    NEXT   J 

9960    NEXT    I 

9965    LPRINT 

9970    RETURN 

10000    REM    CLOSE    ALL  FILES    AND  EXIT    PROGRAM 

10010    PRINT    CHR$(26):PRINT    "DO    YOU   REALLY    WANT   TO    QUIT" 

10  020    INPUT   "ENTER  A  Y  FOR  YES,    AN  N  FOR  NO         ",Z$ 

10030    IF    Z$    =    "Y"    THEN    END    ELSE    RETURN 
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